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Abstract 

We examine the thermal conductivity and bulk viscosity of a one- 
dimensional (ID) chain of particles with cubic- plus-quartic interparticle 
potentials and no on-site potentials. This system is equivalent to the FPU- 
afi system in a subset of its parameter space. We identify three distinct 
frequency regimes which we call the hydrodynamic regime, the perturba- 
tive regime and the collisionless regime. In the lowest frequency regime 
(the hydrodynamic regime) heat is transported ballistically by long wave- 
length sound modes. The model that we use to describe this behaviour 
predicts that as lj ^ the frequency dependent bulk viscosity, ({uj), 
and the frequency dependent thermal conductivity, k{uj), should diverge 
with the same power law dependence on uj. Thus, we can define the bulk 
Prandtl number, Pr^ = fefl("(a;)/(mK(ti;)), where m is the particle mass 
and kB is Boltzmann's constant. This dimensionless ratio should approach 
a constant value as ci; ^ 0. We use mode-coupling theory to predict the 
UJ —> limit of Pr^. Values of Pr^ obtained from simulations are in agree- 
ment with these predictions over a wide range of system parameters. In 
the middle frequency regime, which we call the peHurbative regime, heat 
is transported by sound modes which are damped by four-phonon pro- 
cesses. This regime is characterized by an intermediate-frequency plateau 
in the value of k{uj). We find that the value of k(uj) in this plateau region 
is proportional to where T is the temperature; this is in agreement 
with the expected result of a four-phonon Boltzmann-Peierls equation 
calculation. The Boltzmann-Peierls approach fails, however, to give a 
nonvanishing bulk viscosity for all FPU-a/9 chains. We call the highest 
frequency regime the collisionless regime since at these frequencies the 
observing times are much shorter than the characteristic relaxation times 
of phonons. 

1 Introduction 

The thermal transport properties of one-dimensional (ID) chains of particles 
have been recognized as an enigmatic puzzle for several decades [HI HH] . The 
problem is of more than academic interest because of the speculation that ID 
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systems may have applications for "thermal management" in nanoelectronics 
[32] and it has been suggested that the unusual properties of ID systems might 
make a "thermal transistor" possible [7] . By comparison, momentum transport 
in ID chains has hardly been studied. However, recently there has been work 
[T71[5T] suggesting that thermal transport in ID systems can be best understood 
by considering its coupling to momentum transport. In |17| a theory is developed 
which allows the low frequency part of the heat current power spectrum to be 
predicted from the higher frequency parts of the heat current and momentum 
current power spectra. The theory is shown to provide a good prediction for 
the special case of a system with the ratio of specific heats 7 = cp/cy — 1. 
In such a system the heat current power spectrum can be predicted from the 
momentum current power spectrum alone. One of the purposes of the present 
paper is to demonstrate that the theory developed in [17] also provides good 
predictions for the more general case of 7 ^ 1. 

Fourier's law of heat conduction is Jg — "KVr(r, t) where is the macro- 
scopic heat flux density, ^(r, t) is the local temperature and k, is the thermal 
conductivity. Newton's law of bulk viscous dissipation is Jp = — C(V • v(r, t))fi, 
where Jp is the macroscopic normal momentum current density (or stress) across 
a surface with normal direction n, C, is the bulk viscosity and v(r, is the local 
macroscopic velocity field. We ignore shear viscosity since it is irrelevant for 
ID systems. The sense in which most ID systems do not obey Fourier's law 
is that K fails to converge to a finite macroscopic value [T^]. Rather, k is seen 
to go as 7V^, where TV is the number of particles in the chain and p is some 
positive power. The value of p is a matter of great interest, with different values 
being reported for different systems [HI [191 [29]. Some attempts have been made 
[TT] [2T] to theoretically predict the value of p. While there are arguments in 
favour of universal behaviour [TT] [THl HI] , the existence of universal behaviour 
is not supported by the bulk of simulation results that have been reported so 
far. The consensus understanding of thermal conductivity in ID systems can 
roughly be summarized as follows: the criteria for a system to have a finite 
conductivity are not known, for systems with infinite conductivity the depen- 
dence of the conductivity on system size is not understood and it is not known 
whether this dependence is universal in any sense. 

Momentum transport in ID systems is even more poorly understood, partly 
because it has rarely been discussed in the literature. It has been predicted 
using mode coupling [21 [10] that the momentum current correlation function has 
a long time tail that goes as t~^^'^ in 1D0. This is equivalent to a momentum 
current power spectrum going as for small w or a bulk viscosity that goes 

as N^/'^ . Given that the same calculations predict a i^^/^ time tail for the heat 
current correlation function, which is not supported by any simulation results, 
this prediction should be viewed with scepticism. 

The failure of the mode-coupling calculation in ID is not unexpected. A key 
assumption is that there is a clear distinction between fast and slow relaxation 

^However, as was seen in 1171 . unlike the case for the energy current correlation function, 
the amplitude of the t~l/^ tail in the momentum current correlation function can vanish in 
some systems. 
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processes and that for properties on long time-scales the effects of the fast 
processes can be incorporated into phenomenological hydrodynamic parameters. 
This is roughly correct in 3D where phase space is dominated by large momenta 
and the multiplicity of the large momentum fast processes dominate. In ID the 
low momentum parts of phase space are much more significant and there are 
important relaxation processes on all time scales. A separation into fast and 
slow is still qualitatively correct and perhaps even semi-quantitatively correct 
provided one generalizes the hydrodynamic parameters from fixed constants to 
time-scale dependent ones. 

This picture was tested in [17] for the limiting situation in which mode- 
coupling theory predicted that, because of a vanishing thermodynamic ampli- 
tude, there are no slow momentum current relaxation processes. Here we test 
this picture under more general circumstances. If energy and momentum current 
correlations are dominated by the same fast relaxation processes then mode- 
coupling predicts that these two correlation functions are related by calculable 
thermodynamic quantities. This is a refutable proposition. 

We present simulation results for the FPU-q;/3 model in which the interaction 
has been tuned from nearly harmonic through highly anharmonic and asym- 
metric to ultimately the highly anharmonic but symmetric pure quartic model 
studied in [T7]. The thermodynamic amplitudes predicted by mode-coupling 
theory that we test show a very non-trivial and non-monotonic variation and 
our simulations track this variation exceptionally well in the highly anharmonic 
and asymmetric regime and are not inconsistent at the two extremes. As the 
harmonic limit is approached the comparison with simulation becomes difficult 
because all relaxation processes slow down. At the other extreme, as the pure 
quartic model is approached the slow relaxation part of the momentum current 
correlations becomes too small to observe. 

Because our simulation results confirm the predictions of mode-coupling in 
the limited sense described above, they also lend credence to the "mode cas- 
cade" toy model we proposed in [T7]. However, it is also clear that a complete 
validation of the toy model by simulation will be very difficult because of the 
enormous range of time scales that will have to be tracked. 

The outline for the rest of the paper is as follows. In the remainder of this 
introduction we provide some general background necessary for understanding 
the details of our simulations and comparisons with theory. This is followed in 
section 2 by some specific mode-coupling background. Sections 3 and 4 are a 
description of our model and the thermodynamic calculations we have had to 
perform. Section 5 gives the numerical results of our simulations and is followed 
by a concluding section 6 which is a more extensive discussion than is given 
above. Appendix A gives details on the symplectic algorithm that we are using 
in our molecular dynamics simulations. Appendices B and C discuss issues not 
directly related to mode-coupling but for which our simulations have provided 
insight. We emphasize throughout the paper that it is important to distinguish 
different frequency regimes by the processes that are important, but in practice 
this is not always easy to do as there is, as yet, little guidance from theory. 
Appendix B is a comparison of Boltzmann-Peierls phonon scattering predictions 
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for a weak coupling regime of the FPU-/? model. The theoretical calculations 
are limited to the relaxation-time approximation so that, while the agreement 
with simulation is not perfect, it does suggest that Boltzmann-Peierls can be a 
valid description in a well defined parameter range for the thermal conductivity; 
however, as discussed in Appendix B, the Boltzmann-Peierls approach gives a 
vanishing result for the bulk viscosity. Appendix C discusses the even simpler 
situation of no phonon scattering and is verified, for the momentum current 
correlations, to be reliable for the uppermost frequency range of the FPU-/? 
model. 

The thermal conductivity, /c, and the bulk viscosity, C, are formally related 
to the generalized or wave- vector and frequency dependent transport coefficients 
K{k,u)) and C{k,oj) by [13] 



K = lim k(Lu) = lim lim K{k,Lu) (la) 

C = lim C(^^) = lini lim (^(fc,a;) . (lb) 
Here and in the following we use the notation 

lim^(fc)EEi. (2) 

The frequency dependent k{Lu) and C(^) can be written as Green-Kubo rela- 
tions in terms of the corresponding equilibrium heat current correlation function 
(HCCF), C^it), and momentum current correlation function (MCCF), C(^{t), 
namely 



k(uj) = fdt'e'^''C^{t') (3a) 

t—^OO 2 J — t 

Cicu) = \hnU'dt'e^-''Cdt'), (3b) 

where ks is Boltzmann's constant and /3 = l/ksT is the inverse equilibrium 
temperature. In terms of the corresponding currents j^(t), where fi — n or (, 
we have 

C,{t)^ lun^^(^Sj^{t)Sj,iO)) , (4) 

where L is the system length, (• • •) denotes a canonical average and 5j^(t) is 
the A: — > limit of the deviation of j^ji,{k, t) from its equilibrium value. A further 
remark on the definition of i^{t) is warranted. Following jfi it) is the zero 
A; limit of the current density j^{k,t). Under this definition the heat current 
density, jK{t), is equivalent to what is often referred to as the "total heat flux", 
JK.{t) = J2iLijK{'lijt) {qi is the position of the z**^ particle), in studies where 
Green-Kubo relations are used to calculate thermal conductivities [12]. 
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The HCCF, C^, in the Green-Kubo relation above can freely be exchanged 
with the energy current correlation function, ECCF, [13] defined analogously 
to Ck. but with the heat flux density, j„, replaced with the energy flux density, 
jg. Numerically, it is more convenient to calculate the ECCF and this is what 
we do in the simulations reported in this paper. In practice, we work with the 
Fourier transformed versions or power spectra of the correlation functions, 

~ /"OO 

C^,{uj) = / dte^^[iuot)C^,{t). (5) 

J —oo 

We refer to these as the momentum current power spectrum (MCPS) and the 
energy current power spectrum (ECPS). 

As shown below, in our system the particle spacing is arbitrary. This, as 
well as computational convenience, makes it useful to work in a particle count- 
ing scheme (Lagrange picture) rather than a spatial coordinates scheme (Euler 
picture). In the Euler picture (spatial coordinate scheme) hydrodynamic den- 
sities and currents are expressed as integrals over coordinates containing delta 
functions of the form ^(r — q^) where r is the spatial coordinate being integrated 
over and is the coordinate of the particle labeled i. In contrast, the Lagrange 
picture allows these quantities to be expressed simply as sums over all particles. 
In the Lagrange picture our spatial Fourier transform over a ID chain is defined 
as 

1 ^ 

^j^(fc,i) = ^E'5jrMiK'% (6) 

where k = {-{N - 1)tt/N, . . . , -2tt/N, 0, 27r/iV, . . . , tt} and 5f~^''^ is the devi- 
ation from mean current between particles s and s— 1. The prefactor of iV^^/^ 
is to keep our current independent of system size as an aid to comparison be- 
tween runs. With this definition, {5i^{t)5j^{0)) is not proportional to N and we 
should revise (jl]) by the replacement 1/L — N/L = l/£, where I is the system 
length per particle. However, because t is arbitrary, as will be shown below, it 
is preferable to also define particle based currents (i.e. to work in the Lagrange 
picture instead of the Euler picture). Thus, we use throughout 

C^{t) EE lim (slm^iO)) , (7) 

instead of (|4|). 

To sec how to write a current we examine the energy current. Consider 

the Hamiltonian of the chain 

N 



1=1 



2m 2 ^ 



(8) 



where ^i^j is the potential energy due to the interaction between particles i and 
J. We are currently only interested in the case of nearest neighbour interacions 
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so the sum over potential energies becomes J^iLi'f'i,^-^- Let us divide this 
up between particles so that E = J2i=i -^i where Ei is the energy of the i^^ 
particle. We can choose to divide the energy up in several ways, leading to 
several different definitions of Ei. We choose 

E^ = ^ + l + 0^+1,^] • (9) 

2m 2 

We insist that 4>i,i^i depends only on the difference, q — Qi — Qi-i, between 
the positions qi, qi^i so that d4>i^i-i/ dqi = —d4>i^i-i/ dqi-i. Taking the time 
derivative of Ei we obtain 

d = + 1 [{q, - g,_i)</>-.,_i + - qi)(j}[+i^,] , 

at m at 2 ' 

where 4>'i j{q) = 94'i,j/dq and qi = dqi/dt = Vi. From Hamilton's equations, 

Vi = -</>M-i + so that we get 

Er = \ [{v^+l + - {V^ + I'.- 1 ■ (10) 

Taking the convention that a flow to the right (positive direction) is positive we 
can rewrite (fTU)) as 

E.^~J^"'+]T"\ (11) 
where j\ is the energy flow (to the right) from particle i — 1 to particle i 

.2+1/2 

and je is the energy flow from particle i to particle i + 1. Comparing this 
with (fTII|) we can identify 

if'/' - + v,+x)(j)[+^^, = i(i;, + v,+i)n+i^, . (12) 

where r^+i^^ is the local stress (in ID simply equal to the force of one particle 
on the other) between particles i + 1 and i. A form of the energy current which 
often appears in the literature (e.g. [27j) is 

Jl = ~\W^+l^^+<l^'^,^-l]^^■ (13) 

This definition comes about if one divides up the energy between bonds instead 
of dividing it up between particles (i.e. © is replaced by Ei = [pf ) /4m + 

(/)i+i^i). It is easy to verify that this leads to the same definition of the total 
energy current and so it is equivalent. A similar derivation to the one for (|12p 
leads to a Lagrange picture expression for the momentum current, 

3t"' = -c^,+,, = n+i,. (14) 

Both and HH) are defined in terms of divergences. Hence, they are 
arbitrary up to an additive constant which, if chosen to give zero mean current, 
makes and (5j^ identical. For the momentum current a logical choice for this 
additive constant is the pressure. 
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2 The Bulk Prandtl Number 



It has been known since pioneering studies in the 1960s and 1970s [2 [3] that 
the current correlation functions ([4]) generally have long-time tails that decay 
as some power of time. Mode-coupling theory was developed to explain this 
phenomenon. It invokes a physical picture in which each transport mode in 
the system is coupled to every other transport mode. There are many different 
formalisms of mode coupling theory using many different sets of assumptions. 
A good summary of the most widely used formalisms can be found in [26] . One 
of the most rigorous mode coupling theories is that of Ernst, et. al. p, [TUl[TT| . 
In this theory the main starting assumptions are that the system of interest is 
in local equibrium and that local equilibrium is established and maintained by 
processes which are fast whereas the variations of the hydrodynamic currents 
and the couplings between them occur via processes which are slow. The precise 
meanings of "fast" and "slow" are neither defined as an input to the theory, nor 
are the meanings provided by the theory. Presumably there must be some 
minimum separation in time scales between the fast and slow time scales in 
order for the theory to be valid. 

Of interest to the present study is the prediction, developed in [10] , of the 
leading order terms of the momentum current correlation function and the heat 
current correlation function. These are predicted for a system of general dimen- 
sionality, d. For ID the prediction of TO] can be easily shown [16j for long times 
to be 



where Dt = mn/ pcp is the thermal diffusivity, = (7 — 1)Dt + Di is the 
sound damping constant, k is the thermal conductivity, Di = (477/3 + C)/p 
the longitudinal diffusivity, 77 is the shear viscosity (which vanishes for ID) and 
^ is the bulk viscosity. Here m is the mass per particle, p is the mass density, 
7 = cp/cv is the ratio of the specific heats, cp is the specific heat capacity per 
particle at constant pressure and cy is the specific heat capacity per particle 
at constant volume. The remaining constants in (|15|) are the thermodynamic 
quantities 




(15b) 



(15a) 
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where c is the adiabatic speed of sound, s is the entropy per particle and ap = 
{l/L){dL/dT)p is the thermal expansion coefficient. 

The Green-Kubo integrands in (jlSp have several important features that are 
worth noting. First of all, their power law dependence on t is of great interest; 

the predicted t~^/^ behaviour (equivalent to power spectrum behaviour Ci{u!) ~ 
w"^/^ at small uj) is not in agreement with more recent theoretical [T71 [T51 ?IT\ 
or numerical [T71 [TO] results. Thus, it is tempting to dismiss the prediction 
as entirely incorrect for ID. However, as was shown in .l7j, in ID a careful 
examination of for the case of 7 = 1 yields a prediction that any C'^(t) ~ 
t~^^^ term is zero (raising the possibility of a finite bulk viscosity) and C,^{t) ^ 
and simulations in that same paper confirm this prediction for a chain 
with quartic interparticle potentials. Much more importantly, we showed in [17j 
that by interpreting Ts as a frequency dependent phenomenological parameter 

rs(a;) obtained from simulations of Cq{w) at relatively high frequency, the C'k(w) 
could be predicted and agreed with simulations at low frequency. Thus, it is 
possible that (jisp is only "minimally wrong" ' by which we mean that (fT5a|) 
and (jl5bp are approximately correct if they are interpreted as relations between 
time dependent Tg and Dt for short times and C(; and (7^ at long times. We 
do not specify what constitutes short or long; rather pS]) is to be viewed as 
a pair of renormalization group equations in which every (7^ and Ck. at long 
times becomes the input for the Tg and Dt at the now short times for the 
next iteration. The toy model introduced in the appendix of was meant 

to illustrate this feature. If Ci^{u!) and Ck(u;) have the frequency dependence 
for some range of small uj then over a well defined range of much smaller 
Lo they will vary as lu^'^ with q = (1 — p)/(2 — p). This cascade will eventually 
terminate at the fixed point q = p = p* = (3 — \/5)/2 ~ 0.382. Note that 
(|15|) is a special case of this renormalization group flow; an initial constant Tg 

and Dx implying p = give rise to q = 1/2 which is the dependence and 

Ck oc or equivalently the long time tails cx t^^l"^ in ([TS]). If we now set 

p = 1/2 then the next frequency power becomes q = 1/3 and this is the universal 
behaviour predicted in However, we do not agree with the authors of [3T] 
that there are no further renormalizations. In our view the u:~^l'^ dependence 
of the current power spectra also applies only over a limited frequency interval 
and will, at even lower frequency, be replaced by an uj~'^/^ dependence and 
ultimately hy lo~^ . 

The simulations in [17] tested a very small part of this picture as they relied 

only on the vanishing of and Mhh for 7 = 1 and that = / 0^ . The 

latter simply expresses the fact that heat is carried by weakly damped sound 
modes and does not require the full blown machinery of mode-coupling theory. 
In the present paper we describe much more significant tests. In particular, one 
of the implications of our renormalization group interpretation of (|15p is that, 
at long enough times, T g and Dt will have the same frequency dependence. In 

this case the ratio limt^oo Cc^(t) / Cii{t) = limj^^o Ci;{uj)/Ck.{ui) should approach 
a constant at sufficiently low frequencies. We speculate, given the assumptions 
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in [ini [H] that "sufficiently low frequencies" means on time scales longer than 
those of the processes which establish and maintain local thermal equilibrium. 

This ratio of transport coefficients, which we expect to be constant in the 
thermodynamic and long time limit, is similar to the Prandtl number defined 
as Pr = v/Dt, where v = rj/p is the kinematic shear viscosity. However, our 
ratio is frequency dependent and involves the bulk viscosity. As a convenient 
dimensionless ratio we define the bulk Prandtl number 



Pre ^ lim PrcH ^ lim -^4^ = Itn 



From p5|) and (fTH]) we obtain 



Pre 



f3 



M4 



2Dt) 



1/2 



HH 



which is an implicit equation for Pr^ since 



^ = (7-l) + ^ = (7-l) + |^Prc, 
Ut rnn Kb 



(17) 



(18) 



(19) 



As a first approximation we observe that cp(^/{mK) is normally quite small for 
FPU systems of the type studied here. Assuming this we can write 



approx _ 



Ml 



7-1 



1/2 



M 



(20) 



We can find the exact Pr^ explicitly by noting that (|T8l) and ([T9|) combined are 
a quadratic in Pr^ which can be solved to give 



Pre = + + mJ^M, + ^^Mh + ^ , (21) 

where we must take the positive solution of the quadratic for consistency with 
([20]) in the appropriate limit and where for compactness we have written 

M, = (22a) 
Mh = -^Mhh- (22b) 

It must be stressed that this result should only be correct in a low frequency 
range, below the frequencies of processes which establish and maintain local 
thermodynamic equilibrium. We will see in section [3l over the range of pa- 
rameters where this frequency regime is accessible to simulation, that the above 
provides a good prediction of the bulk Prandtl number in the system of interest. 
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3 The Cubic-Plus-Quartic Chain 



We consider a ID chain system with nearest neighbour forces governed by a 
cubic-plus-quartic interparticle potential. The Hamihonian is 



w(p,q) = E 



N r 2 



T, y ^{(li - ~ a) + —{q,- q,^i- a) 

2m 6 4 



(23) 



where pi = mdqi/dt is the momentum of the i^^ particle, qi is the position of 
the i^^ particle, m is the mass of each particle, a is the equilibrium interparticle 
spacing which is a fixed, arbitrary, system parameter and a and B are force 
constants. We use periodic boundary conditions so that qo = qN — L. This is 
the familiar FPU-a/? system with the harmonic force constant set to zero. We 
use B for the coefficient of the quartic term in the potential so as not to confuse 
it with the inverse temperature, j3 = (ksT)'^. We also choose to work with 
interparticle distance as coordinates, rather than absolute particle position, so 
we define Xi = qi — qi^i — a. 

There are a number of statistical mechanical quantities of interest to us in 
this work and these are most easily calculated in an isobaric ensemble. The 
statistical weight of a configuration is proportional to the Boltzmann factor 
modified by a pressure term and L = qj^ — qo = Na + Xi is allowed to take 
on all possible values so that the partition function takes the form [28] 



exp(-/3G) = -j^ J dpdqexp[-(3in + PL)] 



1 N 



dp / dXF{[3,P) 



(24) 



where dpdq denotes the phase space volume element for particles i — 1 . . .N ,h\s 
Planck's constant, G — G(/3, P, N) is the Gibbs free energy and P is the pressure 
(same as average force in ID). We note that there is no need for a prefactor of 
in the partition function because each particle is connected only to its 
neighbours so that particle ordering makes the particles distinguishable. The 
N statistically independent factors, F{(3, P) are 



F(/3, P) = exp 



2m 3 4 



PX 



(25) 



We will restrict our chain simulations to zero pressure but the formal expression 
(j25p that includes P ^ is convenient for deriving various thermodynamic 
quantities. By defining characteristic length and time scales and using these to 
express the dynamical variables X and p in dimensionless form we can eliminate 
much of the apparent parameter dependence in (|25p . We choose length and time 
units 
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4 = (/3B)-^/^ to^(^^y\ (26) 

and write 

X = eox, p='^v, (27) 

where now x and v are dimensionless. Our Boltzmann factor in these new 
dimensionless variables is 

Fif3, P) = exp (^-^ - a*^ - ^ - P*x^ , (28) 

which shows that the parameter space of the model is two-dimensional with two 
dimensionless parameters defined by 



P* = pP£o (29) 

= ivT- (30) 

We are currently interested in the zero pressure system. This leaves us with a 
parameter space of interest that is one-dimensional and we can study the whole 
parameter space by fixing any two of a, /9, and B and varying the third. Because 
we have previously studied the pure quartic case [T^ it is convenient for us to 
fix (3 — 1, B = 1 and vary a. As can be seen in (PU)) . varying a is equivalent 
to varying the temperature, T. To summarize, our units are defined by m = 1, 
B = 1 and = 1. 

The a* = case is the pure quartic chain studied in [17j while in the limit 
of a* — !■ oo the system approaches the harmonic chain. This can be seen by 
expanding the potential, V, about its minimum. With reference to ()23|) we 
define 

This has its minimum at Xcq — ~a*. Expanding around the minimum using 
X — Xcq + Sx we can write the dimensionless, scaled potential as 

f3V{6x) = const. + ^Sx^ - '^^^^ + ■ (32) 

If l|32p is used in the Boltzmann factor the harmonic term limits the magnitude 
of Sx to 0{l/a*). Thus, the cubic and quartic terms will be 0(l/a*^) and 
0(l/a*^) respectively and, hence, negligible in the a* oo limit. Between the 
extremes of a* = and a* — *■ oo the potential is both highly anharmonic and 
asymmetric and this is precisely what we need for a significant test of the Ernst 
formulae (fT5|) in a case distinct from that examined in [l^. The fact that by 
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varying a* we can also explore the crossover from strong to weak anharmonicity 
is an added bonus. Given the above discussion, removing the constraint of 
zero harmonic force constant in (|23p would probably not yield much additional 
information. 

We are now ready to discuss the various quantities involved in calculating 
Pre;: c, cp, 7, Dt, Ts, cup, dap/dT, dcp/dT, dc/dp. We will make use of the 
constant pressure partition function (j24[) which reduces to 



exp (-/3G) = 



27rm\ 



N/2 



^-N/3Pa 



/oo 
dX exp [-/? {V{X) + PX)] 
-OO 



(33) 



with V{X) given by pip . In deriving the thermodynamic quantities of interest 
we repeatedly encounter averages of the form 



_ r°° dXX''exp\-f](V(X) + PX)] 

X„(/3, P) = 4"5„ = (X"> = ^^§5 f ^ r . \ ^ > (34) 

^ ' " " \ / J^^dXexp[-(3{V{X) + PX)] ' ^ ^ 

which is the n'^ moment of X . These moments satisfy the recursion relation 

(n + 1)X„ = (3 [aX„+3 + BX„+4 + PX^+i] • (35) 

Hence, besides the trivial Xq — xq = 1, we need only evaluate Xi and X2 as 
numerical integrals. Furthermore, the temperature and pressure derivatives of 
'Xn follow directly from ([34]). We get 



-^Xn = -13 (X„+i - X„Xi) , (36a) 

+ ^{XnXi - Xn+i) + P(XnXi - X„+i) . (36b) 

The relations p6ap and (|36b[l form the basis for our subsequent derivation 
of the thermodynamic quantities. For future reference we list below a number of 
the most important thermodynamic expressions characterizing the cubic-plus- 
quartic chain. As a useful starting point, we obtain an equation of state via 

{L) = — ^N{a + X,{(3,P)). (37) 
and hence our equilibrium spacing between neighbouring particles is 

i= ^ = a + Xi^a + ioxi . (38) 
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The fact that l|55|) depends on the arbitrary spacing, a, impHes that £ can- 
not appear in any expression for a thermodynamic quantity except as a triv- 
ial multipUer. The average energy per particle is found from the enthalpy, 

(E) + P{L) = d{l3G)/dl3, and is 



TV 



-2^3 



1 



3 " 4 

We obtain the thermal expansion coefficient by differentiation of ([38|) 



(39) 
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1 1 /df 



1 1 



-P{X2 - Xi 



T£ 
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— (X4 - X3X1) + -{X5 - X4X1) + P*{X2 - Xl^) 



(40) 



Similarly, the isothermal compressibility is 
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XT 
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The specific heat capacity per particle at constant pressure is 



(41) 



cp 



Cp_ 

N 



N 'dp 



m+p{L)] 



kf 



1 a* _ 2\ ot* ,_ , 1 / 



X/^ 



2P*a* P* 

^ {X4 - X3X1) + — {X5~ X4X1) 



(42) 



and for the specific heat at constant volume one can use the identity cy = 
Cp — £Tap^/xT with the results from (|^ and (|¥T|) . The specific heat ratio 7 = 
Cp/cv, combined with (|4ip gives the adiabatic speed of sound c = {1/ pXsY^"^ , 
where Xs = Xt/7 is the adiabatic compressibility, as 



mXT 



X2 



(43) 
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As noted above, the sound speed varies linearly with i as expected; this is its 
only dependence on the arbitrary spacing, a. In fact the arbitrariness is absent 
in the Lagrangian picture which we adopt for our simulations; here the natural 
unit is c/£ and the speed of sound is viewed as being measured in terms of 
particle number per unit time (i.e. it is a "hopping" frequency for a disturbance 
to move from one particle to the next). 

At zero pressure for both limits a* = and a* — > cx), the specific heat 
ratio 7=1. Numerical evaluation of the formulae above shows the maximum 
7 ~ 1.544 at a* ~ 2.418. This value of a* serves as a reasonable central value 
for our chain simulations. 

4 Evaluation of the Bulk Prandtl Number 

We wish to predict Pr^ as a function of a* . In order to do this we need the 
thermodynamic expressions from the previous section. We must also be able to 
evaluate the Ernst current amplitudes (|16p . Let us now reduce the expressions 
of the current amplitudes to a form which is useful for numerical evaluation for 

the cubic-plus-quartic model. is already in a convenient form. We note 

that ^ is clearly proportional to the mean lattice spacing, ^, as is expected 

of a quantity proportional to a current. But this same linear dependence is not 
evident in and Mhh- We rewrite them as 



We now express mnn and m_| in terms of the quantities c/£ and £ap which do 

not depend on the arbitrary spacing parameter, a. Doing this, and noting the 
structures of ((44|) . which involves expressions of the form (l/A){dA/dT) and 
{l/A){dA/dp) we are able to express them as 



The specific heat relation cp—cy — iTap'^/xr or {j—l)/{£ap) = apjT/ {xrcp) 
allows us to eliminate the potentially vanishing denominators in (|45p . so that 
these relations remain valid in the limits where 7 — )■ 1 and ap — > 0. We replace 




(44b) 



(44a) 




(45b) 



(45a) 



([45bll by 
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"f^HH apjT /dln{£ap/cp)\ 



I XTCP \ dT J 



(46) 



p 



and note that this can be evaluated by the procedures described in section [3l 
These procedures do not apply directly to the derivative at constant s in (j45ap 
but we can start with the chain rule and write 



(47) 



The coefficient {dP/d£)s of the second term in (gZl) is -l/{£xs) = -7/(^Xt), 
which is expressed entirely in terms of quantities found in section |31 The coef- 
ficient of the first term can be rewritten 



ap^T 

XTCp 



(48) 



in which we have used the Maxwell relation {dS/dt)T = {dP/dT)i. Thus, in 
summary, we can calculate by 



iapT fdlii{T£/c)\ (d\n{£/c 



CP \ dT J p ' \ dP ''^^^ 

which is expressed entirely in terms of the quantities found in section [3l The 
derivative calculations are very messy and algebraic packages such as Maple 

are very useful to ensure correctness of the final results. With and tuhh 

in hand we can now proceed to calculate Pr^ for any value of a* using (PT|) . 
We thus obtain the theoretical curve shown in Fig. [T] As discussed earlier, 
we are restricting our attention to the region in parameter space surrounding 
a* ~ 2.418. At this value of a* the value of 7 is maximized. We see in Fig. [1] 
that Ptq has a pronounced local minimum at a* ~ 2.418. 

The values of Pri^ obtained from simulations, described in the next section, 
can be seen in Fig. [1] as well. The simulation results were obtained by varying 
a with kpT = 1, B = 1, P = m = 1. As will be discussed, the simulations 

indicate that Cq{uj) and Ck.{ijj) go as the same power of w as (jJ — > and so we 
find a well defined value of Prc^ at each value of a* . As can be seen in Fig. [TJ 
the agreement between the simulations and the theory is good. This supports 
the physical picture that we have been proposing. 

5 Numerical Results 

We carry out molecular dynamics simulations of the system described by (|23p in 
periodic boundary conditions for various values of the dimensionless parameter 
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Figure 1: The value of the bulk Prandtl number, Pr^, as a function of the 
dimensionless parameter a* . We have fixed i? = 1, /3 = 1, to = 1, P = 0, and 
varied a* by varying the cubic coefficient a. The simulation points for a* = 
1.5, and 3.2 are shown even though, at the lowest frequencies examined in those 
simulations, the hydrodynamic regime had not yet been reached. 
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Figure 2: Current power spectra for a representative sample of values of a*, 
a) a* = (pure quartic), b) a* = 2.0, c) a* = 2.418, d) a* = 2.8. The 

theoretically predicted curves for were obtained by the method described in 
[TT] . The approach of the observed energy power spectrum to the theoretical 
curve occurs at lower frequency as the value of a* is increased. 
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a*. We choose N — 2^^ and typically run for 2^^ time units, outputing 4 times 
per time unit. We set m — 1, ksT — 1, B — 1 and a = 0. We vary a* by 
changing a. This, combined with our Hamiltonian (j23p . defines our units of 
length, time and energy. In particular, distance and time are measured in units 
defined by l|26p. Each run is initialized by randomly generating the positions 
and momenta of the particles according to an isobaric ensemble. With the 
algorithms that we use we can achieve good accuracy using 8 time steps per 
time unit. We output the sums over all particles of the momentum current 
and energy current. Then, using the methods outlined in [T7], we calculate 
the momentum and energy current power spectra. In this work we have used 
eighth order, sixth order and fourth order symplectic algorithms. The eighth 
order algorithm was used for the a* ~ 0, a* ~ 2.0 and a* — 2.418 cases. 
It is a refinement of the algorithms presented in [31]. This algorithm is the 
same as the one used in [17]. It was realized later in this work that relaxing 
the symmetry used by Yoshida to obtain this eighth order algorithm allows the 
development of sixth order algorithms which are, nevertheless, more accurate 
than the eighth order one in terms of the error in energy. Further, although 
the energy error was increased significantly by using a fourth order symplectic 
algorithm there was still no secular change in energy. As a result, there is no 
appreciable loss in accuracy when a fourth order symplectic algorithm is used 
to find the total momentum and energy currents as we do in this work. Hence, 
the later simulations in this work (a* > 1.8) were carried out using fourth or 
sixth order symplectic integrators. These were checked against small numbers 
of runs using the eighth order routine to verify that the accuracy was preserved. 
The coefficients of our fourth and sixth order integrators appear in Appendix 
A. 

From our simulations we thus obtain Ci^{uj) and Ce{uj) (recall that in the 
Green-Kubo equation for k{uj) we can freely replace Ck. with Ce). We also use 
the theory developed in section II of [l7] to produce a theoretical prediction of 

Ce{i^). Varying a* we see that there appear to be distinct regimes in which the 
current power spectra behave in different ways. Plots representative of these 
regimes are presented in Fig. [2] For reasons discussed in detail in [17] the a* = 

case is special with C(^{uj) const and with (^^(a;) — > as w ^ 0. For all 

nonzero values of a* , Cq{uj) — + cxd as ~> 0. We note that the observed C^{u!) 
always approaches the theoretical prediction, but that this approach happens 
at progressively lower frequencies as a* is increased. 

Let us examine the different frequency regimes present in the transport in 
these systems. In the highest frequency regime, defined by lot >> 1 where r 
is a typical mode lifetime, damping plays no role . Accordingly we call this 
the collisionless regime. This frequency regime will be examined in Appendix 
C. At very small frequencies the hydrodynamic mode coupling model presented 
in [17j accurately predicts the energy current power spectrum. This suggests 
ballistic transport of heat via low frequency sound modes which damp via a 
sound damping coefhcient which is frequency dependent. Further evidence that 
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Figure 3: The frequency dependent bulk Prandtl number as obtained by simu- 
lation of the cubic-plus-quartic chain with a* = 2.0. The horizontal dashed line 
is our asymptotic estimate and is given as a point in Figure [TJ 

this is occuring in this regime is that, for a* ^ 0, and go as the same 
power of Lo in this regime. This is consistent with the picture presented in our 
toy model presented in [17j . This can be interpreted as evidence that local 
equilibrium is established and maintained by fast processes and gives rise to the 
mode coupling behaviour, as is assumed in lOj. We call this frequency regime 
the hydrodynamic regime. 

As we increase a* we see that an intermediate frequency regime becomes 

established. This regime has a C'c(aj) ^ u)""^ part which gives way to a constant 

Cf:{oj) plateau. This is characteristic of the existence of a single relaxation time 
in this regime. We note that as a* increases the interparticle potential becomes 
more harmonic in the vicinity of its minimum. Thus we might expect that, in 
some frequency regime, the behaviour of the system might look more and more 
like a harmonic lattice with small anharmonic perturbations as we increase a*. 
This resembles the well known Boltzmann-Peierls picture of heat transport in 
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which heat is carried at the sound velocity c by phonons which scatter over 
some mean free path. In the absence of defects the scattering is entirely due 
to phonon-phonon interactions. It is known that, at least at first order, for ID 
systems this damping cannot be due to the cubic term (three phonon processes) 
in the potential [TH j so that the damping must be due to the quartic part 
of the potential (four phonon processes). It is worth noting, however, that this 
argument is only valid if the cubic and quartic coefficients are small. If the 
coefficients are as given in ([5^ then the second order contribution due to three 
phonon processes can be comparable in size to the first order contribution due to 
four phonon processes. In any case, given our speculation that the Boltzmann- 
Peierls picture is correct in the intermediate frequency regime for this system we 
refer to this intermediate frequency regime as the perturbative regime. We give 
further evidence below that this picture is correct, though it should be stressed 
what we may be seeing is superposition or interference of second order cubic 
with first order quartic effects. We discuss this regime still further in Appendix 
B where we work with a pure quartic model to avoid the difficulties of the second 
order cubic contributions. A puzzle is that we see only a single relaxation time 
in the perturbative regime. 

It is simple to estimate Pri^{uj) from our simulation results by calculating 

Cq/ {m(3Cf) for each frequency. A typical example is shown in Fig. [3l At the 
lowest frequencies examined in the simulation Pr^(a;) seems to have converged 
to a constant value. This is indicative of the hydrodynamic regime being reached 
at these frequencies. For each value of a* we take the average value of Pr^{u!) 
in the frequency range for which it appears to have converged and use this 
value as an estimate of Pr^. The resulting values of Pr,^ for the values of a* 
examined are seen in Fig. [TJ For values of a* both higher and lower than the 
range examined, the frequency regime in which PrQ{uj) converges to a finite 
value is pushed to frequencies too low to be easily accessible. The simulation 
points are shown for a* — 1.5 and a* ~ 3.2 even though the simulations did not 
probe low enough frequencies to see the hydrodynamic regime at those values of 
a* . Thus the disagreement seen in Figure [1] for these two points should not be 
construed as a failure of the theory to predict Pr^. Similar comments apply to 
simulations carried out for a* — 3.6, 4.0, 4.2, 4.4 and 4.6. Indeed, via methods 
from [17] we can predict the frequency at which we should see the crossover 
to the hydrodynamic regime and for a* = 3.2 the crossover should occur at 
uj/2tt ~ 2~^^, so the lack of agreement seen for a* > 3.2 is expected. For 
this reason we have not shown results for a* > 3.2 in Figure [TJ In the range 
from a* = 1.8 to a* = 2.8, where the hydrodynamic regime was accessible, 
the agreement between simulation and the prediction from ()21[) is good. The 
standard errors in the simulation values are too small to be shown. However, 
significant systematic errors are expected to be present since we are taking the 
average value of PrQ{uj) which is going asymptotically to a constant value. In 
particular, the value of PrQ{uj) may not have been very close to its asymptotic 
value for the case of a* — 2.418. The Pr,^ vs. uj curve showed a discernable 
slope at the lowest frequencies examined for this value of a* . 
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Figure 4: The current power spectra for the cubic-plus-quartic chain with a* = 
3.2. The hydrodynamic regime is at lower frequencies than are probed in this 
simulation. The perturbative regime is very evident with the low frequencies 
following a Lorentzian shape. A Lortentzian fit to the energy current power 

spectrum is shown. The spectrum is well fit by Cf* = C'e(0)/(1 + w^/wq) with 

(7,(0) = 2.03 X 103, LOO = 6.91 x IQ--*. 
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Let us examine the perturbative regime where the behaviour resembles damp- 
ing with a single relaxation time. This regime covers a wider frequency range as 
we increase a* beyond about 3.0. An example of this is shown in Fig.[5]in which 

a Lorentzian fit is shown to Cc{lo). We must stress that this fit is not the — > 

limit of Ce{uj), since the hydrodynamic regime is at lower frequencies than those 
shown in Fig.d] If this behaviour applied as w ^ then we would see a finite k. 
Hence, if we ignore the behaviour in the lower frequency hydrodynamic regime, 
we may define an effective thermal conductivity in the single relaxation time 
regime which is just the zero frequency value that we obtain by a Lorentian fit 
to the heat current power spectrum in this regime. This effective, perturbative 

regime, thermal conductivity is Kq'''^*"'^'^ = [Cf^iuj — 0)]/2. We can ask how 
this n depends on the temperature. Noting that we can think of varying a* as 
equivalent to varying T we can define an effective temperature, T*. From pop . 
if B = 1 and a = 1, a* - T-^/^. Thus, we define 

T* = (a*)-* . (50) 

A plot of Kq'^'^'"'^'' vs. T* , is shown in Fig. [5] and wc see that for low temper- 
ature Kq'''*"'^'^ goes as (T*)"^. This is consistent with the well known result of a 
four phonon Boltzmann equation [20l[24]. We see this as strong evidence that, 
for approximately harmonic lattices, there is a perturbative regime at frequen- 
cies higher than the hydrodynamic regime. Whether this perturbative regime is 
well described by the four-phonon Boltzmann equation in every detail remains 
to be checked; a limited comparison is given in Appendix B. 

In Appendix B we also indicate why a Boltzmann-Peierls equation approach 
fails to predict the chain bulk viscosity. 

6 Conclusions 

Several key points are worth pointing out immediately. 

I. For the cubic- plus-quartic system, in general, there are at least three dis- 
tinct frequency regimes. The lowest frequency regime, which we call the 
hydrodynamic regime, is characterized by ballistic transport of heat via 
long wavelength sound waves. The next regime, which we call the per- 
turbative regime resembles a harmonic chain with a quartic perturbation 
in some respects. In this regime the transport of heat can be viewed as 
being damped by four phonon processes; the damping appears to be gov- 
erned by a single relaxation time in the cubic-plus-quartic system but by 
a range of relaxation times in the absence of a cubic term as described in 
Appendix B. At the highest frequencies observing times are too short for 
any appreciable phonon scattering. We call this the collisionless regime 
and we examine it further in Appendix C. A striking aspect of these results 
is that for nearly harmonic systems the perturbative regime and the hy- 
drodynamic regime are clearly separated. For more anharmonic systems 
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Figure 5: The value of the perturbative regime plateau value of the thermal 

conductivity, Kq*^'^*"'^'', as a function of the effective temperature T* for low 
temperature (large a*). We have fixed B = 1, P = 1, m = 1, and varied a* 
by varying the cubic coefficient a. The dotted line is for reference and shows 
(T*)-2 behaviour. 
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the perturbative regime is not visible because the hydrodynamic regime 
is estabhshed at higher frequencies. 

2. The theory developed in |il7] predicts the thermal conductivity well in the 
hydrodynamic regime. This theory is based on an assumption of heat be- 
ing ballistically carried by sound waves. Thus, it is reasonable to conclude 
that the momentum transport and heat transport are strongly coupled. 
This, combined with the fact that the power spectra have the same power 
law behaviour, is evidence that a theory like the toy model from the ap- 
pendix of |17| describes the transport of both quantities. The current 
paper shows that, using the results of [IQJ, we can extend the theory in 
[T7j to predict the bulk viscosity as well. 

3. The regime that we have called the perturbative regime is characterized 
by a K ~ behaviour. This is indicative of damping of modes by 
4-phonon scattering (or likely a combination of first order effects from 4- 
phonon scattering and second order effects from 3-phonon scattering) as 
one would calculate using the four-phonon Boltzmann equation [201 121] ■ 
However, this identification is speculative. 

4. For a system like the one examined here with 7^1 the bulk viscosity 
is infinite. At sufficiently low frequencies the momentum current power 
spectrum has the same power law behaviour as the heat current power 
spectrum. 

5. In the hydrodynamic regime the mode-coupling theory of Ernst et. al. 
[SI dni dlj fails to correctly predict the correct power law divergence of the 
transport coefficients. However, it does predict the correct ratio between 
the thermal conductivity and the bulk viscosity (the bulk Prantdl num- 
ber). This theory assumes that local thermal equilibrium is maintained 
by fast processes whereas the hydrodynamic transport of heat and mo- 
mentum is carried out by slower processes. The success of the theory in 
predicting the bulk Prandtl number is evidence that the assumptions of 
this theory are correct for a ID anharmonic chain with 7 7^ 1. Further, 
our results allow us to quantify what is meant by "fast" and "slow" in 
this system. Specifically, the mode cascade of our model from [T7j shows 
that on any time scale within the hydrodynamic regime the slow relax- 
ation processes on that time scale are produced by much faster processes. 
Thus, the meaning of "slow" is simply the time scale on which we are 
observing the relaxation process. The meaning of "fast" is the frequencies 
for which Ffe/fc'^ ~ w as is discussed following equation (10) in [T7]. Thus, 
the meanings of "fast" and "slow" vary with the time scale that we are 
examining. 

We have, in this paper and in [TT, described heat in ID systems as being 
ballistically transported by sound waves which are weakly damped. Purely bal- 
listic transport would have no damping at all and would exhibit Ce{t) t~^. 
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Ballistic transport of heat resembles second sound in which heat is transported 
at constant speed (the speed of ordinary sound in ID); an important difference 
is that second sound exists only over very restricted temperature and frequency 
windows 0] whereas ballistic transport is not subject to these restrictions. It 
might seem to be more appropriate to refer to the transport seen here as su- 
perdiffusive transport {Ce{t) ~ t~P, p < 1). However, we feel that the term 
"superdiffusive transport" implies a process similar to diffusion (a random walk 
where an initial distribution spreads with a speed which decreases with time) 
whereas what we are describing is much closer to ballistic transport in which 
a "packet" of heat is carried at constant speed over some long distance by a 
sound wave. Further, if the toy model of [T7] is correct then this produces a 
mode cascade with an infinite series of different power law behaviours of the 
transport coefficients. This is quite distinct from the usual picture of superdif- 
fusion. Thus, we propose the term "damped ballistic transport" for this type of 
transport. 

The existence of the perturbative and hydrodynamic regimes is highly sig- 
nificant when interpreting our results and the results of others. It is possible for 
the plateau in the perturbative regime to look as if k(w) is converging to a finite 
value at w ^ 0. However, if lower frequencies (larger systems) were examined 
the hydrodynamic regime would be seen and k(Lu) would be seen to diverge. 
Further, the extremely low frequency at which the hydrodynamic regime as- 
serts itself in the system examined in this paper is below the lowest frequencies 
examined throughout much of the literature (for a summary of relevant results 
see [19]). As always, extreme caution must be exercised in claiming that one is 
seeing the thermodynamic limit. 

Another interesting feature of these results is that it seems to be valid to 
examine frequencies lower than the fundamental frequency of the system. That 
is, we can observe processes on time scales longer than £/c, the time for a long 
wavelength disturbance to travel around the system one time. Indeed, we can 
change the size of the system while holding our simulation times constant and see 
no change in the lowest frequency behaviour of the system even for frequencies 
below the fundamental frequency. This is another consequence of the mode 
cascade. For any system size we can observe relaxation processes to time-scales 
that the fundamental frequency dominates, and these slowest observed processes 
are much slower than the fundamental mode oscillations. The lowest frequency 
transport modes that we are observing propagate around the system many times 
over the course of the simulation, but since their damping is driven by much 
faster processes it does not matter that they are circulating around a short chain 
rather than passing once through a very long chain. 

There are several areas of further work which are needed to complete the 
above picture. While the 4-phonon Boltzmann equation calculation has been 
done for the steady state case [35] , no such calculation has been carried out for 
the equilibrium state. This calculation should be carried out and the overall 
amplitude of the damping in the perturbative regime needs to be calculated so 
that the four-phonon Boltzmann equation can be used to quantitatively predict 
the heat current power spectrum in the perturbative regime. Also, we do not 



25 



have a satisfactory explanation for why, in the perturbative regime, a single 
relaxation time dominates. This question might be answered by an equilibrium 
calculation. A summary of the results of the nonequilibrium calculation in [22j 
and analysis of how these results can be interpreted in the context of this study 
is presented in Appendix B. Also discussed in Appendix B is the failure of 
the Boltzmann-Peierls approach for the calculation of the chain bulk viscosity. 
Clearly a completely new approach is required for this property. 

To show that the physical picture proposed in our toy model is correct, 
simulations would need to be carried out which probe systems to low enough 
frequencies to observe the onset of the next power law behaviour in the power 
spectra. We predict that for the cubic-plus-quartic system at a* = 2.0 the 
"kinks" in the power spectra corresponding to the onset of the next power law 
behaviour in the series should occur somewhere below uj ~ 2~^^. For now, such 
low frequencies are beyond the practical limits of our simulations. However, 
some other system might demonstrate the transitions between power laws at 
more accessible frequencies. 

A Symplectic integrators 

There is by now a very large literature on improvements to the Yoshida in- 
tegrators and a particularly extensive list of alternatives is given by [53]. We 
have not tested or used any of the versions that require calculations of the force 
gradient in addition to the force. Instead we have opted for those schemes in 
which improvements are obtained solely by relaxing the requirement that the 
integrator contain the minimum possible number of steps. This can reduce, in 
some cases dramatically, the size of the intermediate steps in the integrator and 
will typically improve both stability and accuracy at a small increase in running 
time. 

We have tested a limited number of integrators and give below two inte- 
grators, one 4th order and one 6th order, that we have found nearly optimal 
and used in production runs. Both are "position" integrators which means the 
first move is a position update. This is followed by a momentum update, then 
position, in the alternating sequence 



Xnew = w(l)At, Pnew = P + fw{2)/S.t , 

m 
P 

Xnew = X -\ w{3)At, Pnew ^ P+ fw{4:)At . . . , (51) 

m 

where / is the force. The coefficients are as follows. 

4th order integration coefficients 
w{l)=w{9) = 5/3/{3 + V39), w;(2) = zz;(8) - 3/4 , 
w(3) =w(7) = -2/3/(6 -f%/39), w(4) u;(6) = -1/4 , ^ ' 

w(5) =23/3/(4-1- VM) 
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(53) 



6th order integration coefficients 
0.055, 

0.f52f2924f8f 98012208708832, 
0.2381129197090666122795175, 
0.3450759351638895426864652, 
0.5086956937118671097820404, 
-0.0432317055841977505550037, 
-0.4259353129298358880967277, 
0.0460265286005069869976553, 
0.2482533990178043320703396 

The 4th order coefficients are very close to the optimal equation (62) coeffi- 
cients in j23j but have the advantage, as analytical expressions, of being easily 
programmed for arbitrary precision. The 6th order coefficients are based on 
the equation (83) "momentum" integrator of |23j but have an added position 
step at the beginning which has been adjusted to further improve the integra- 
tor. It is worth remarking that finding new solutions to the required non-linear 
constraint equations defining a symplcctic integrator is a very difficult search 
problem but that modifying an existing one by small parameter increments is 
easy using Newton-Raphson iteration. 
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B Phonons and the Boltzmann-Peierls Equation 

Our ultimate goal is to test the accuracy of the Boltzmann-Peierls approach to 
thermal transport in ID chains in the frequency regime between the hydrody- 
namic regime discussed in the main sections of this paper and the collisionless 
regime discussed in Appendix C. However, at the present time the theory is 
not well enough developed to make definitive tests possible. In this appendix 
we report on the more modest achievement of a comparison of simulation with 
Boltzmann-Peierls in the relaxation time approximation. Even this has only 
become possible because of the recent study ^22, by one of us (BN) of the 
Boltzmann-Peierls equation for a chain of weakly anharmonic oscillators in a 
thermal gradient. An intermediate result in that study was an analytical for- 
mula that can be used to predict the wave-vector dependent relaxation rate of 
phonons in the relaxation time approximation. This in turn allows us to predict 
the energy current power spectrum which we compare to simulations based on 
the model used in |22|. The results are good enough to unambiguously ver- 
ify that the low frequency structure we see in the simulations is indeed due to 
phonon relaxation and that the Boltzmann-Peierls picture is qualitatively and 
even semi-quantitatively correct. 

The introductory part of our discussion relies heavily on results in [22 . To 
avoid excessive duplication we will use the same conventions as were adopted 
in that paper. The relevant model is the FPU-/? model, namely equal mass 
particles described by the classical Hamiltonian similar to ([23|) 
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E 



' (54) 



where the sum extends over the N particles in the chain with Xi+ n = Xi while 
Xi is the deviation of a nearest neighbour pair separation from the equilibrium 
spacing a and is defined as below (j23p . In the weakly anharmonic limit we can 
use the harmonic part of 7i to define the normal modes of the system. These 
are the phonons labelled by the N independent wave- vectors fc„ = 2-Kn/N and 
the phonon frequency 



sin (— ) 



UJZB^^J^ (55) 
V m 



has its maximum ojzb the Brillouin zone boundary while the phonon group 
velocity Uk = dujk/dk vanishes there. The term in (|54p gives rise to phonon 
scattering and is treated in the Boltzmann-Peierls approximation in which the 
phonon occupation number is treated as a local density and its rate of change 
by transport, drik/dt + Ukdrik/dx, is equated to the rate of change by collisions, 
Rki^ki)- In [22] only the transport term Ukdrik/dx was relevant as the system 
was a steady state, non-equilibrium oscillator chain in a thermal gradient. Here, 
for describing the relaxation of fluctuations from equilibrium in an otherwise 
spatially homogeneous chain, only duk/dt is relevant. As in [22j we write the 
deviation from equilibrium in terms of a new function defined by 

Suk = Uk- ""^Uk = ""^nkCuk + l)gk , (56) 

where "^^n^ is the Bose factor 1/ {exp(fiujk / ksT) — 1). The transport term in the 
Boltzmann-Peierls equation can then be written 

dnk eq ,eq , -, n %fe /rrrN 

nk{ uk + l)^, (57) 
and this is to be equated to the net collision rate which to linear order in g is 



Rkig) = -TT-^^Tri / dki / dk2UJk^kiUJk2^^k3 

ion K2 J J 

X {9k - 9ki - 9k2 + 9k3) ""^rik ''W^ 

xCWi + l)CW2 + l)'5(wfc -Wfei -iOk2 +t^fc3)' (58) 

where the fci, ^2 integrations are understood to be over an interval of 27r and 
ks = ki + k2 — k mod 2tt. The rate (|58|) is based on Fermi's golden rule and its 
derivation is standard textbook material that we need not repeat here. 
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The general solution to "''nkCnk + l)dgk/dt = Rk{g) can in principle be 
given in terms of the solutions to the associated eigenvalue equation but to our 
knowledge these have never been obtained and we do not attempt such solution 
here. Instead we make what is known as the relaxation time approximation 
which is to set the deviations gfc^, gk2, and 17^3 in the integrand in ([55)) to 
zero. Then Rk{g) becomes gk times an integral over equilibrium distributions 
which we can choose to write as the wave- vector dependent function — "^^rife i^'' 
Uk + ^)/Tk- That is, with this approximation, the solutions to the Boltzmann- 
Peierls equation are 



5fc oc exp (-i/Tfe) , (59) 



with 



Tk IGtt" K2^ 

eq eq / eq , a \ / e.q , i \ 

X V nk^i ^rifci + 1)( ^nk2 + 1) 

X 5 (wfe - Wfci - + ujk3 ) ■ (60) 

Note that we are only interested in the classical limit in which case we can set 
^''uk ^'^'^rik + 1 ~ kBT/fnok and the relaxation rate equation (|60p becomes 



1 = risin^fJjX;,, (61a) 



Tk \2 



9 fK^ksT 
16^"^^ [-1^2 



Kk = dki dk2UJzBS{uJk - t^ki - i^k2 + ^ka) , (61c) 



where we have used also (jSSp to express uJk in terms of wave-vector k. 
The dimensionless integral Kk was evaluated in [52] with the result 



Kk = 2 0^ {z-i/'^B(l/3,l/3)F(l/3,l/3;2/3,;z) 

-2^/^5(2/3, 2/3)F(2/3, 2/3; 4/3; z)| , (62a) 

where B{x,y) — r(x)r(y)/r(a; + y) is the beta function and the F — 2F1 are 
hypergeometric functions. 

The relaxation time approximation deserves some comment. Were we to 
make the corresponding approximation of setting gk^ = gk2 — Sfea = in the 
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thermal gradient problem in [22], we would find the amplitude of the leading 
divergence in the phonon distribution to be wrong by factor four. However 
in that problem we know the distributions are driven by the thermal gradient 
and therefore the gk at different k are strongly correlated. In the equilibrium 
fluctuation problem we are considering here, it is quite reasonable to suppose 
that the gfc. (fc^ ^ k) in (|58p are only weakly correlated with 17^ and thus 
on average nearly zero even when is not. This makes the relaxation time 
approximation much less severe here; admittedly we cannot conclude anything 
quantitative. 

The independent exponential relaxation in time of each normal mode as 
given in l|59p implies the energy current power spectrum is a mode sum of 
the corresponding frequency domain lorenztians (x (l/Tfe)/(cj^ + (1/rfc)^). The 
explicit formula is 



C,{u) = 2{kBTf / —Uk \ , , , (63) 



2^"'= (c.2 + (l/rfe)2) 

which also expresses the fact that energy in mode k is transported at the group 
velocity — {1 / 2)uj z b cos{k / 2) . The temperature dependent prefactor in (I63p 

can be deduced from the sum-rule requirement that J{duj/2T:)C^{uj) is the equal 
time energy current-current average < Sj^ > which in the harmonic limit is 
{kBTYijJzB^ /9>. Given the explicit formulae l|61a[) for the relaxation time we 

can write Cc{u}) in the scaling form 

d,{u:) = {kBTf'^d\uj/^), (64) 
where C^{uj) is the dimensionless spectrum 



C~M.j-.„'y-y^'y-^^--^,. ,65, 

The integral in (|65p must be done numerically but then, as given by (I64p . can 
be applied universally with only an amplitude and frequency rescaling. The 

spectrum C^iuj) is even in lo, has the normalization 

dujC^{u) = I , (66) 

and is characterized by simple power-law behaviour in the two limits of low and 
high frequency. These are, for positive 



C^{uj) « tj-2/5Q 2749172... ^ (t^^o) (67a) 
w w^^i 1488360, (t^ ^ 00). (67b) 

The Lo^"^/^ low frequency behaviour arises because of the small k singularity in 
the relaxation time integral given in (|62ap . This in turn, as discussed in j22| . 
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Figure 6: Energy current power spectrum deviation from best fit in the form 
5 = log2[C'e(w)S™"i'^"™/C'e(w)Ei'^inj] vs. log2(w/27r). The deviations have 
been shifted for clarity. The relaxation rate correction factors for the curves 
from low to high A are Rr = 0.607, 0.893, 1.114, 1.317, 1.455, 1.547, 1.704, 
1.803, 1.847, 1.910, 1.960. The fit spectrum ^ is characterized by the appro- 
priately scaled asymptotes (|67ap and (j67bp : their intersection defines a corner 
frequency Wcornor = 2.4443i?ri^ = i?T-2.4443(9/87r)A~''/2 that is shown as the 
heavy vertical lines. 
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is a consequence of the fact that the phonon dispersion curve tUk at small k is 
linear with a cubic correction. Thus the u;~^/^ divergence observed here is quite 
general and not at all special to the nearest neighbour model ([51)1. On the other 
hand, it may require the absence of odd terms in the potential. The evidence for 
this is the frequency independent regime seen in the cubic-plus-quartic model 
spectrum in Fig. 3] 

For purposes of comparing the above theoretical results to numerical simula- 
tion we first note that we can rescale lengths and times such that the Boltzmann 
factor exp(— H/fcsT) is reduced to a one parameter family. The specific scaling 
we have chosen turns Ti{p,X) in into 



V^^n[v,x) = ^1+a\x^ + \x\ (68a) 

and is equivalent to setting m — K4 ~ ksT = 1, K2 — A. Comparison with ([32]) 
shows that, except for the cubic term, the FPU-/? model here is the (x-shifted) 
FPU-a/3 model of the text with A — (a*)^. The case A = is the pure quartic 
model discussed in [T7] whereas weak anharmonicity requires A 00. The 
largest A in the simulations described below is A = 20. 

A number of thermodynamic properties of the model can be given explicitly 
in terms of known functions and are useful as numerical checks. First note that 
because the potential terms in (|68ap are even in x, the specific heat ratio 7 = 1. 
We can then get the adiabatic sound speed, c, trivially from the isothermal 
compressibility and find 



^ - iF-T- (69a) 

where if 3/4 and K1/4 are Bessel functions and cq^ = A is the result one would 
have in the absence of the term in (j68ap . Similarly, for the equal time energy 
current-current average we get 

^R,o 4' > / < Sj,^ >o= I [3i?A' - 1] (70) 

where < Sj^ >o~ A/2 is the result with no x"^ term in (|68ap and equals 

^ {duj / 2'k)C ^{lo) in the relaxation time approximation. Thus any deviation of 
the ratio Rjo in (|70p from unity indicates a failure of the power spectrum sum- 
rule. This failure vanishes as ^ ^ 00 but is about 3% at ^ = 10 and rises 
to 11% at A = 5 and 27% A ~ i. To some extent the failure is spurious 
and just a consequence of our having chosen to determine the amplitude in (|63p 
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Figure 7: Individual Rr estimates from Fig [6] vs. 1/A. The solid box near the 
top of the graph is the data from runs on chains of length 2^^; all other data 
points are for N — 2^^ . The data is well summarized for A > 3 by the smooth 
curve which is Rr = (2.056 - O.IO/A^ - A.m/A^)/{1 + 20.47/^^) ^nd provides 
an estimate of the weak coupling limit A ^ oo from our runs at finite A. We 
use units with Ki{= B) = 1^ fi = \ and m = 1. 
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by evaluating averages in the harmonic Hmit. A perfectly reasonable expedi- 
ent which we adopt is to redefine C^{lu) = 2AC^{uj /VL) as given in to 

C^{u)) =^ Rjo2AC^{uj and this ensures the sum-rule is exactly satisfied. 
We also recognize that the relaxation time approximation for the mode decay 
rate, namely l/r^ = fl sm^{k/2)Kk from (|61ap . is unlikely to be exactly correct. 
The very simplest correction one can make here is a single multiplicative con- 
stant for all 1/rfe and hence we make the replacement O —^ Rt^. In conclusion, 
we compare 

4,.) . („) 

to the simulation results with R^- as adjustable and designated as the relaxation 
rate correction factor. 

Comparison with simulation is shown in Fig. [6] in the form of a logarithmic 

deviation, namely log2[C'e(w)S'"™i'^"°"/C£(u;)'^'i"l^] vs. log2(t^/27r). Values for 
A range from 3 to 20 and it is only for the smallest A that any significant 
systematic deviation in the low frequency phonon region can be detected. Each 
fit in Fig. [6] has yielded a relaxation rate correction factor and these are shown 
in Fig.[7|as Rt vs. 1/A. Any variation with A indicates that there are processes 
contributing to the scattering beyond that given by the Fermi golden rule result 
((60|) and the data in Fig. [7] is consistent with these perturbative corrections 
scaling asymptotically as l/A^. This is expected since, as given in (j68bp . 1/A 
is proportional to ^/Kl in the original Hamiltonian and thus the corrections 
vary linearly with perturbation, K4. Most significantly, the simulation results 
show that there is a limiting value for Rr as A —> 00 and, thus, are consistent 
with the scaling expected from the Boltzmann-Peierls analysis. The fact that 
this limiting value for Rr is not unity but Rr w 2.0 — 2.1 means our analysis is 
not yet truly quantitative. We have no way of knowing whether the observed 
rate correction factor is the result of the relaxation time approximation, failure 
of the Boltzmann "Stosszahlansatz" in the context of ID phonons, or some 
combination of the two. Resolution will have to wait until further theoretical 
work is completed. 

The above discussion suggests that we could calculate the corresponding 
chain viscosity from the Boltzmann-Peierls equation. We indicate briefly why 
this approach fails. The physical picture underlying the calculation of viscous 
dissipation in (3D) insulating solids is due to Akhieser |T] - see also [HI [HI 130] ■ 
Imagine one of the low frequency phonons in a lattice system of phonons. What 
is meant by low frequency in this case is lut < 1, where r is the average relax- 
ation time (due to phonon collisions) of all phonons in the system. The given 
phonon can be thought of as slowly modulating the spatial density. Therefore, 
because of anharmonicity in the interparticle interactions, the elastic constants 
and frequencies of all phonons are also modulated by the phonon in question. 
The mode Gruneisen parameter 7^, where for ID [5] 
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7fc = - 



{L) dujk 
ujk d{L) 



a(ln(L)))' 



d{lnuJk) 



(72) 



describes the frequency shift of mode k due to the lattice strain induced by the 
given phonon. In general (3D) the 7^ vary from mode to mode, so that the 
differing frequency shifts put the system temporarily out of thermal equilibrium 
(recall that the thermal mode occupancies are adiabatic invariants and thus 
temporarily constant - i.e. until relaxation occurs - for small w, and that the 
new Uk will depend on the shifted frequencies). If lu is small enough (i.e. lot <C 
1) the system can relax back to equilibrium during the period of modulation, 
211/10, of the lattice and the distortion process is quasi-static (i.e. reversible 
or dissipationless) . For larger w, say lot < 1, the relaxation will not precisely 
follow the distortion, leading to irreversible behaviour with viscous dissipation. 

For frequencies lot < 1 Akhieser and others [U [H [151 EQ] have applied the 
Boltzmann-Peierls equation, with and without the single relaxation time ap- 
proximation, to calculate 3D lattice viscosities (i.e. shear and bulk). If we 
apply this theory to ID chains we find that the bulk viscosity is a sum of con- 
tributions from all modes, and that the contribution of mode k contains the 
factor {Sjk)'^ , where 6jk — Ik ~ (7)7 with (...) denoting an average over all 
k. As we shall see, for our purposes the exact nature of the averaging over k 
need not be specified. The appearance of the quantity 5^k is unsurprising in 
view of the physical picture given above. We now evaluate 5^k for a ID chain 
with Hamiltonian of type ([54]) and ([23]) . with K3, and K4 denoting the 
harmonic, cubic and quartic spring constants, respectively. The value of for 
such a chain is [5] 



Note that 7/c is independent of k, so that S'jk = 0. Thus the Boltzmann-Peierls 
approach fails to predict a viscosity for all such chains (with nearest-neighbour 
interactions), for arbitrary values of K3 and K^. 

C Momentum current in the collisionless regime 

We verify statements made on pages ITSl and [22l of the text that, at least for a 
specific model, the high frequency weak coupling regime is a distinct collision- 
less regime in which phonon scattering plays no role. Explicitly, we show that 
the momentum current correlation function for the FPU-/3 model of Appendix 
B, evaluated as an average over a thermal population of undamped phonons, 
agrees with numerical simulation. Moreover, the comparison with the numer- 
ical simulation shows in what frequency regime other processes dominate the 
momentum current power spectrum. 

The momentum current for the FPU-/3 model of Appendix B is given by 




(73) 
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and it is a straightforward, albeit tedious, textbook exercise to express the co- 
ordinate deviations in terms of phonon creation and annihilation operators 
and evaluate the average {Sj(;{t)Sj^{0)) in an ensemble of non-interacting, Bose 
distributed phonons. The terms proportional to and K2K4 are time inde- 
pendent because the coordinate sums in (fM)) dictate that there are no remaining 
fc 7^ phonons in the final expressions. The term proportional to K4^^ is a three 
phonon average in which the total momentum ki + k2 + ~ 0. The final ex- 
pression for the temporal Fourier transform of the correlation function, in the 
high temperature or classical limit, is 

d^{^) = - (^^^] I dki f dk2 y di±LJi ± W2 ± W3 - tv). (75) 

TT \mLOzB^ J J J ^ 

The integrals in (|75p are understood to be over an interval of 27r and the remain- 
ing sum is meant to indicate eight separate terms corresponding to all possible 
sign combinations in the (5-function. All of these terms can be combined into 
one by the expedient of a factor two, dropping the absolute value signs in the 
definition (j55p of the phonon frequency, and extending the integration interval 
for each momentum to 47r. The result is our coUisionless power spectrum 

= lh^I]SlA-^I{oj/u:zB) (76) 

TT LOZB 

with A given by (|68bp and I{z) the dimensionless spectrum 



I{z) = 1/4 dki dk2S 



ki\ , . fk2\ . fh k2 
sml- +sm - -sm - + - 



(77) 



in which integration intervals of An are understood making I(z) a full period 
integral over a periodic function. 

The evaluation of I{z) will be given below but a number of important features 
can be understood without detailed calculation. A key observation is that (|77p 
is completely analogous to a density of states expression commonly found in 
solid state physics and as such has singular points which are the well known van 
Hove singularities. The three sine function sum in the (5-function has several 
symmetry related local maxima which are also the absolute maxima. One of 
these is at fci = ^2 = 47r/3 and this implies the spectrum approaches a constant 
as z approaches 3-\/3/2 from below and vanishes identically for z > 3\/3/2. Since 
I{z) is an even function of z, there is a corresponding cutoff at z = —3^/3/2. The 
integrand region near the origin, fci = fc2 = 0, is also the source for a van Hove 
singularity and this is rather more unusual and interesting. Because the linear 
terms in the sine sum vanish, the leading behaviour is the cubic proportional to 
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kik2{ki +^2) and simple power counting then shows the spectrum must diverge 
as We can expect this to be a general result, dependent only on the 

fact that phonon dispersion curves vary linearly with cubic corrections at small 
fc, and not specific to the nearest neighbour force constant model being treated 
here. 

The evaluation of (|77p is quite involved with a lot of similarities to what was 
done in [22] to obtain the relaxation time integral Kk that has been reproduced 
here as (|62ap . The first step is the change of variables ki — 2{u+v), k2 = 2{u—v) 
which transforms the sine sum into 2 sin(u) cos(w) — sin(2u) and makes the v 
integration trivial. The remaining u integral is 



-y/4sin^(u) - (sin(2M) + z)^ 

with the u integration over that interval between and tt for which the argu- 
ment of the square root is non-negative. The substitution x = cot(u) puts the 
integrand in (|78p into algebraic form and shows the result can be expressed as an 
elliptic integral. The most convenient result however is obtained by expanding 
the elliptic integral as a series and recognizing the series as the expansion of 



I(z) 



4\ 1/3 



I{z) = - B(l/3,l/3)i^(l/3,l/3;2/3;4z727) 



~t (i) 2/3)^(2/3, 2/3; 4/3; 4^2/27) (79) 

which, surprisingly, except for normalization and a change of argument 42:^/27 
2sin2(fc/2)/27, is the Kk integral ([5^ . The final proof that ^ is the correct 
value of (|78p is completed by showing both forms satisfy the same differential 
equation. All of the above manipulations are complicated and could not have 
been done without computer packages such as Maple. 

The power spectrum (j76p with the explicit result (j79p has been compared 
with the numerical simulations described in Appendix B. The agreement is excel- 
lent for the largest A values and can be improved for smaller A by two renormal- 
izations similar to what was done for the energy current power spectrum in Ap- 
pendix B. First, the harmonic model zone boundary frequency ojzb — 2^/ K2/rn 
is rescaled by the sound velocity ratio c/co using ([69|l . This guarantees that 
oj = u}{k) in the limit /c ^ is given exactly even when, as A is decreased, the 
modes at large k strongly damp and their frequency becomes ambiguous. Sec- 
ond, the spectrum amplitude is rescaled to give the exact frequency sum-rule. 

The integral over the coUisionless Cq{lo) given in (|75|) is 

dwCc(w) = ^4 (80) 

whereas the exact result, as obtained from the equal time thermodynamic av- 
erage 2tt {5 jfj, is 
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Figure 8: Comparison of simulations of the momentum current power spectrum 
with the coUisionless spectrum ([7S)) and ([75]) rescaled as described in the text. 
The horizontal line is the estimated amplitude for the zero frequency limit of 
the spectrum for the pure quartic model taken from [17j and is based on sim- 
ulations on chains longer than those treated here. Arrows on curves for ^4 > 
mark the expected frequency for the longest possible wavelength mode in the 
system. Even harmonics of this fundamental are apparent in the large A, weakly 
anharmonic, simulations; the fundamental itself is not visible except at small 
A. As in Fig. [SI the simulations are run with K4 — f3 — m — 1, K2 — A. 
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dujCi^{uj) — TrkBTK2 



3Rk - 1 



RkA^ 



(81) 



where Rk is the Bessel function ratio (j69bp . The exact expression reduces, in 
the two hmiting cases of large and smaU A to 



\2TTkBTK2 

■nkBTK2 
A 



18 
l2 



12Rt - 



Rr 



0{A) 



where 



Rr 



r(3/4) 



(82) 



(83) 



r(i/4) 

Thus we confirm (j80[) in the hmit A oo and yet obtain a finite resuh for 
A^O. 

Comparisons of the rescaled predictions with simulation are shown in Fig. [51 
Remarkably, even the A = 0, pure quartic, model is in qualitative agreement 
at high frequencies with the coUisionless approximation. At lower frequencies 
we see another contribution to the simulation spectrum which presumably is 
related to the finite lifetime of the phonons. Exactly what this relation might 
be is, however, not obvious since, as shown in Appendix B, the Boltzmann- 
Peierls approach will not yield a non-zero viscosity. Also, the corner frequencies 
at which the spectrum saturates in the momentum current power spectrum are 
much lower than those in the energy current spectrum as seen by comparison 
of Figures [8] and [6l Another striking observation from Fig. [8] is that the low 
frequency saturation value of the momentum current spectrum is independent 
of the magnitude of the anharmonic term in the hamiltonian. Such universality 
begs for an explanation; we do not have one. 
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